# Direct outputs eligible for public release

# GLOBAL SETTINGS --------------------------------------------------------------

options(
    scipen = 999,
    digits = 16,
    max.print = .Machine$integer.max,
    show.error.locations = TRUE,
    warn = 1
)

RNGkind("L'Ecuyer-CMRG")
seed <- 818675309L
set.seed(seed) # setting main seed

# PACKAGES ---------------------------------------------------------------------

library(data.table)

library(devtools)
devtools::install_github("setzler/eventStudy/eventStudy")
library(eventStudy) # for DiD-IV

# PACKAGE SETTINGS -------------------------------------------------------------

# data.table
setDTthreads(threads = 1)
options(datatable.print.class = TRUE, datatable.print.keys = TRUE)
# so that printing the data.table also shows the variable type on top

# BEGIN FILE -------------------------------------------------------------------

# Just in case something masks 'source()', use extra deliberate base::source().
base::source("~/code/0-utility-functions/table-and-text-utils.R", local = TRUE)

# DECLARE DATA FROM OUR AND COMPARISON STUDIES ---------------------------------

#' Annuity factor
#'
#' @param r Interest rate.
#' @param delta Time preference rate.
#' @param t Time relative to K. Lottery win occurs at t = 0.
#' @param K Age, average age of winner in our sample: 46
#' @param T Maximum age, 80.
#'
annuity_factor <- function(r, delta, t, K = 46, T = 80) {
    ((1 + r) / (1 + delta))^t * (delta / (1 + delta)) * (1 - (1 / (1 + delta))^(T - K + 1))^(-1) # nolint
}

#' Lambda in IRS (2001)
#'
#' @param r Interest rate.
#' @param delta Time preference rate.
#' @param t Time relative to K. Lottery win occurs at t = 0.
#' @param K Age, average age of winner in our sample: 46
#' @param T Maximum age, 80.
#'
lambda <- function(r, delta, t, K = 50, T = 80) {
    sum_1 <- 0
    sum_2 <- 0
    for (s in 1:20) sum_1 <- sum_1 + (1 + r)^(-s)
    for (s in 1:(T - K)) sum_2 <- sum_2 + (1 + delta)^(-s)
    ((1 + r) / (1 + delta))^t * (sum_1 / sum_2)
}

comp_dt <- list()

# IRS (2001)
# Respondents to a Survey in Massachusetts (Sample Excluding Large Winners)
# This is the headline estimate mentioned in their abstract, and based on their
# preferred sample of winners with annual pre-tax prize payouts of less than
# \$100,000. The estimate we report is based on their Table 4, Column VIII.
# Following their own suggestion in footnote 20, this estimate must be divided
# by $\lambda = 0.9$ in order adjust for the fact that lottery payouts in their
# study are for 20 years instead of an annuity.
comp_dt[[1]] <-
    data.table(
        study = "IRS (2001)",
        src_pop_1 = "Respondents to a Survey",
        src_pop_2 = "in Massachusetts",
        src_pop_3 = "(Sample Excluding Large Winners)",
        posttax_hh_total = as.numeric(NA),
        posttax_winner_total = as.numeric(NA),
        posttax_winner_total_d015 = as.numeric(NA),
        posttax_winner_net_d015 = as.numeric(NA),
        pretax_hh_total = as.numeric(NA),
        pretax_winner_wage = as.numeric(NA),
        pretax_winner_wage_d10 = -1 * (0.101 / 0.9)
    )

# BFGI (2021)
# Administrative Data on Parents with Young to Early-Adult Children
# Do not report income effects in their paper.
comp_dt[[2]] <-
    data.table(
        study = "BFGI (2021)",
        src_pop_1 = "Administrative Data on",
        src_pop_2 = "Parents with Young to",
        src_pop_3 = "Early-Adult Children",
        posttax_hh_total = as.numeric(NA),
        posttax_winner_total = as.numeric(NA),
        posttax_winner_total_d015 = as.numeric(NA),
        posttax_winner_net_d015 = as.numeric(NA),
        pretax_hh_total = as.numeric(NA),
        pretax_winner_wage = as.numeric(NA),
        pretax_winner_wage_d10 = as.numeric(NA)
    )

# Our estimates
# In Table 6.2, Columns 3, 4, 7 and 8, we report MPE estimates based on the
# IV model introduced in Section 3.2 and defined by equations (3.5) and (3.6).
# We obtain these estimates by 2SLS estimation of the two-equation system with
# the endogenous variable being the  (pre- or post-tax) measure of unearned
# income in a given period, and the outcome variable being a (pre- or post-tax)
# measure of individual or household total labor earnings. To calculate unearned
# income, we annuitize the (pre- or post-tax) lottery winnings as explained in
# Section 4.1, assuming  that both the interest rate $r$ and discount rate
# $\delta$  are 2.5 percent, and that the total lifespan is 80 years.
comp_dt[[3]] <-
    data.table(
        study = "Our estimate",
        src_pop_1 = "Administrative Data on",
        src_pop_2 = "Working-Age Individuals",
        src_pop_3 = "and Households",
        posttax_hh_total = -0.523,
        posttax_winner_total = -0.433,
        posttax_winner_total_d015 = as.numeric(NA),
        posttax_winner_net_d015 = as.numeric(NA),
        pretax_hh_total = -0.332,
        pretax_winner_wage = -0.241,
        pretax_winner_wage_d10 = as.numeric(NA)
    )

# The remaining columns 5,6, and 9 are populated with values that we can infer
# from the available estimates, but adjusted for the fact that the discount rate
# is different than the one we assumed in our annuitization. To adjust our
# estimates, we follow IRS (2001) and CLNÖ (2017) and assume that the
# per-period utility function takes the Stone-Geary form. As explained in the
# appendix, we then obtain adjustment factors.
a01 <- vector(mode = "numeric", length = 80 - 46 + 1)
# annuity factor when r == delta, average age in our data is 46
a11 <- vector(mode = "numeric", length = 80 - 46 + 1)
# annuity factor when r != delta
for (t in 0:(length(a01) - 1)) {
    a01[t + 1] <-
        annuity_factor(r = 0.025, delta = 0.025, t = t, K = 46, T = 80)
}
t <- NULL
rm(t)
for (t in 0:(length(a11) - 1)) {
    a11[t + 1] <-
        annuity_factor(r = 0.025, delta = 0.015, t = t, K = 46, T = 80)
}
t <- NULL
rm(t)
adj1 <- mean(a01[2:6]) / mean(a11[2:6])
# adjustment factor, five years following lottery win
comp_dt[[3]][, posttax_winner_total_d015 := posttax_winner_total * adj1]

a02 <- vector(mode = "numeric", length = 80 - 46 + 1)
# annuity factor when r == delta, average age in our data is 46
a12 <- vector(mode = "numeric", length = 80 - 46 + 1)
# annuity factor when r != delta
for (t in 0:(length(a02) - 1)) {
    a02[t + 1] <-
        annuity_factor(r = 0.025, delta = 0.025, t = t, K = 46, T = 80)
}
t <- NULL
rm(t)
for (t in 0:(length(a12) - 1)) {
    a12[t + 1] <-
        annuity_factor(r = 0.025, delta = 0.015, t = t, K = 46, T = 80)
}
t <- NULL
rm(t)
adj2 <- mean(a02[2:6]) / mean(a12[2:6])
# adjustment factor, five years following lottery win
comp_dt[[3]][, posttax_winner_net_d015 := -0.360 * adj2]

a03 <- vector(mode = "numeric", length = 80 - 46 + 1)
# annuity factor when r == delta, average age in our data is 46
a13 <- vector(mode = "numeric", length = 80 - 46 + 1)
# annuity factor when r != delta
for (t in 0:(length(a03) - 1)) {
    a03[t + 1] <-
        annuity_factor(r = 0.025, delta = 0.025, t = t, K = 46, T = 80)
}
t <- NULL
rm(t)
for (t in 0:(length(a13) - 1)) {
    a13[t + 1] <-
        annuity_factor(r = 0.025, delta = 0.100, t = t, K = 46, T = 80)
}
t <- NULL
rm(t)
adj3 <- mean(a03[2:6]) / mean(a13[2:6])
# adjustment factor, five years following lottery win
comp_dt[[3]][, pretax_winner_wage_d10 := pretax_winner_wage * adj3]

# CLNÖ (2017): The estimates corresponds to their estimate in their Table 6,
# run various specifications (for details see ./matlab/estimate.m)
# requires matlab; if not available, will skip to saved results
tryCatch(
    {
        system('matlab -nodisplay -nosplash -nodesktop -r "run ~/code/4-analyses/comparison/matlab/main.m; exit"') # nolint
    },
    error = function(cond) {
        message("Errored; probably MATLAB issue. Skipping to stored results")
    },
    warning = function(cond) {
        message("Warning. Probably don't have MATLAB installed. Skipping to stored results") # nolint
    }
)
matlab_dt <- fread("~/code/4-analyses/comparison/matlab/output/table.csv")

comp_dt[[4]] <-
    data.table(
        study = "CLNÖ (2017)",
        src_pop_1 = "Administrative Data on",
        src_pop_2 = "Working-Age Individuals",
        src_pop_3 = "and Households",
        posttax_hh_total = matlab_dt[model == "household-pretax-2.5%", mpe], # pre-tax here refers to the wage measure, not the lottery win #nolint
        posttax_winner_total = matlab_dt[model == "winner-pretax-2.5%", mpe],
        posttax_winner_total_d015 = matlab_dt[model == "winner-pretax-1.5%", mpe], # nolint
        posttax_winner_net_d015 = matlab_dt[model == "clno", mpe], # CLNO's original specification #nolint
        pretax_hh_total = as.numeric(NA),
        pretax_winner_wage = as.numeric(NA),
        pretax_winner_wage_d10 = as.numeric(NA)
    )

# PSO (2018)
# Do not report income effects in their paper.
comp_dt[[5]] <-
    data.table(
        study = "PSO (2018)",
        src_pop_1 = "Administrative Data on",
        src_pop_2 = "Working-Age Individuals",
        src_pop_3 = "and Households",
        posttax_hh_total = as.numeric(NA),
        posttax_winner_total = as.numeric(NA),
        posttax_winner_total_d015 = as.numeric(NA),
        posttax_winner_net_d015 = as.numeric(NA),
        pretax_hh_total = as.numeric(NA),
        pretax_winner_wage = as.numeric(NA),
        pretax_winner_wage_d10 = as.numeric(NA)
    )

comp_dt <- rbindlist(comp_dt, use.names = TRUE)

# Housekeeping
a01 <- NULL
a11 <- NULL
adj1 <- NULL
a02 <- NULL
a12 <- NULL
adj2 <- NULL
a03 <- NULL
a13 <- NULL
adj3 <- NULL
matlab_dt <- NULL
rm(a01, a11, adj1, a02, a12, adj2, a03, a13, adj3, matlab_dt)

# MAKE TABLE -------------------------------------------------------------------

filedir <- "~/paper/tables/table-6.2.tex"
file.create(filedir)

# Start Table 6.2
start_latex_table(
    column_structure = "lccccccccc",
    file_dir = filedir,
    append = TRUE
)
add_latex_table_rule(rule_type = "toprule", file_dir = filedir, append = TRUE)
add_latex_table_rule(rule_type = "midrule", file_dir = filedir, append = TRUE)

panel_a_heading_input1a <-
    "&  & \\multicolumn{4}{c}{Post-Tax Unearned Income} & & \\multicolumn{3}{c}{Pre-Tax Unearned Income}" # nolint
write_latex_table_row(
    panel_a_heading_input1a,
    file_dir = filedir,
    append = TRUE
)
add_latex_cmidrule(
    point_start_col = "l",
    point_end_col = "r",
    start_col = 3,
    end_col = 6,
    file_dir = filedir,
    append = TRUE
)
add_latex_cmidrule(
    point_start_col = "l",
    point_end_col = "r",
    start_col = 8,
    end_col = 10,
    file_dir = filedir,
    append = TRUE
)

panel_a_heading_input1b <-
    "& Primary Data Source and & Household & Winner & Winner & Winner & & Household & Winner & Winner" # nolint
write_latex_table_row(
    panel_a_heading_input1b,
    file_dir = filedir,
    append = TRUE
)

panel_a_heading_input1b <-
    "Study & Population of Study & Total Labor Earnings & Total Labor Earnings & Total Labor Earnings & Net Labor Income & & Total Labor Earnings & Wage Earnings & Wage Earnings" # nolint
write_latex_table_row(
    panel_a_heading_input1b,
    file_dir = filedir,
    append = TRUE
)

panel_a_heading_input1c <-
    " & & discount rate: 2.5\\% & discount rate: 2.5\\% & discount rate: 1.5\\%  & discount rate: 1.5\\% & & discount rate: 2.5\\% & discount rate: 2.5\\% & discount rate: 10\\%" # nolint
write_latex_table_row(
    panel_a_heading_input1c,
    file_dir = filedir,
    append = TRUE
)

panel_a_heading_input1d <-
    "(1) & (2) & (3) & (4) & (5) & (6) & & (7) & (8) & (9)" # nolint
write_latex_table_row(
    panel_a_heading_input1d,
    file_dir = filedir,
    append = TRUE
)

add_latex_table_rule(rule_type = "midrule", file_dir = filedir, append = TRUE)

# Add body of Table 6.2
# - introduce convenience functions given a similar structure throughout

Table62Row1 <- function(dt, text) {
    dt_temp <- copy(dt)
    dt_input <-
        sprintf(
            "\\multirow{3}{*}{%s} & %s & \\multirow{3}{*}{%s} & \\multirow{3}{*}{%s} & \\multirow{3}{*}{%s} & \\multirow{3}{*}{%s} & & \\multirow{3}{*}{%s} & \\multirow{3}{*}{%s} & \\multirow{3}{*}{%s}", # nolint
            dt_temp[study == text]$study,
            dt_temp[study == text]$src_pop_1,
            fmt_fixed_decimal(
                input = dt_temp[study == text]$posttax_hh_total,
                sub_for_missing = "NA",
                should_be_rounded = TRUE,
                round_digits = 3,
                decimal_digits = 3,
                scale_factor = NA,
                add_comma = TRUE
            ),
            fmt_fixed_decimal(
                input = dt_temp[study == text]$posttax_winner_total,
                sub_for_missing = "NA",
                should_be_rounded = TRUE,
                round_digits = 3,
                decimal_digits = 3,
                scale_factor = NA,
                add_comma = TRUE
            ),
            fmt_fixed_decimal(
                input = dt_temp[study == text]$posttax_winner_total_d015,
                sub_for_missing = "NA",
                should_be_rounded = TRUE,
                round_digits = 3,
                decimal_digits = 3,
                scale_factor = NA,
                add_comma = TRUE
            ),
            fmt_fixed_decimal(
                input = dt_temp[study == text]$posttax_winner_net_d015,
                sub_for_missing = "NA",
                should_be_rounded = TRUE,
                round_digits = 3,
                decimal_digits = 3,
                scale_factor = NA,
                add_comma = TRUE
            ),
            fmt_fixed_decimal(
                input = dt_temp[study == text]$pretax_hh_total,
                sub_for_missing = "NA",
                should_be_rounded = TRUE,
                round_digits = 3,
                decimal_digits = 3,
                scale_factor = NA,
                add_comma = TRUE
            ),
            fmt_fixed_decimal(
                input = dt_temp[study == text]$pretax_winner_wage,
                sub_for_missing = "NA",
                should_be_rounded = TRUE,
                round_digits = 3,
                decimal_digits = 3,
                scale_factor = NA,
                add_comma = TRUE
            ),
            fmt_fixed_decimal(
                input = dt_temp[study == text]$pretax_winner_wage_d10,
                sub_for_missing = "NA",
                should_be_rounded = TRUE,
                round_digits = 3,
                decimal_digits = 3,
                scale_factor = NA,
                add_comma = TRUE
            )
        )
    write_latex_table_row(
        dt_input,
        file_dir = filedir,
        append = TRUE
    )
    invisible(dt_input)
}

Table62Row2 <- function(dt, text) {
    dt_temp <- copy(dt)
    dt_input <-
        sprintf(
            " & %s &  &  &  & & & &  & ", # nolint
            dt_temp[study == text]$src_pop_2
        )
    write_latex_table_row(
        dt_input,
        file_dir = filedir,
        append = TRUE
    )
    invisible(dt_input)
}

Table62Row3 <- function(dt, text) {
    dt_temp <- copy(dt)
    dt_input <-
        sprintf(
            " & %s &  &  &  & & & &  & ", # nolint
            dt_temp[study == text]$src_pop_3
        )
    write_latex_table_row(
        dt_input,
        file_dir = filedir,
        append = TRUE
    )
    invisible(dt_input)
}

write_latex_table_row(
    "\\emph{Panel A: U.S. Studies} & & & & & & & & ",
    file_dir = filedir,
    append = TRUE
)
write_latex_table_row(
    " & & & & & & & & & ",
    file_dir = filedir,
    short_row_height = TRUE,
    append = TRUE
)

Table62Row1(dt = comp_dt, text = "IRS (2001)")
Table62Row2(dt = comp_dt, text = "IRS (2001)")
Table62Row3(dt = comp_dt, text = "IRS (2001)")
write_latex_table_row(
    " & & & & & & & & & ",
    file_dir = filedir,
    short_row_height = TRUE,
    append = TRUE
)

Table62Row1(dt = comp_dt, text = "BFGI (2021)")
Table62Row2(dt = comp_dt, text = "BFGI (2021)")
Table62Row3(dt = comp_dt, text = "BFGI (2021)")
write_latex_table_row(
    " & & & & & & & & & ",
    file_dir = filedir,
    short_row_height = TRUE,
    append = TRUE
)

Table62Row1(dt = comp_dt, text = "Our estimate")
Table62Row2(dt = comp_dt, text = "Our estimate")
Table62Row3(dt = comp_dt, text = "Our estimate")
write_latex_table_row(
    " & & & & & & & & & ",
    file_dir = filedir,
    short_row_height = TRUE,
    append = TRUE
)

add_latex_table_rule(rule_type = "midrule", file_dir = filedir, append = TRUE)

write_latex_table_row(
    "\\emph{Panel B: European Studies} & & & & & & & & ",
    file_dir = filedir,
    append = TRUE
)
write_latex_table_row(
    " & & & & & & & & & ",
    file_dir = filedir,
    short_row_height = TRUE,
    append = TRUE
)

Table62Row1(dt = comp_dt, text = "CLNÖ (2017)")
Table62Row2(dt = comp_dt, text = "CLNÖ (2017)")
Table62Row3(dt = comp_dt, text = "CLNÖ (2017)")
write_latex_table_row(
    " & & & & & & & & & ",
    file_dir = filedir,
    short_row_height = TRUE,
    append = TRUE
)

Table62Row1(dt = comp_dt, text = "PSO (2018)")
Table62Row2(dt = comp_dt, text = "PSO (2018)")
Table62Row3(dt = comp_dt, text = "PSO (2018)")
write_latex_table_row(
    " & & & & & & & & & ",
    file_dir = filedir,
    short_row_height = TRUE,
    append = TRUE
)

# End Table 6.2
add_latex_table_rule(
    rule_type = "bottomrule",
    file_dir = filedir,
    append = TRUE
)
end_latex_table(file_dir = filedir, append = TRUE)
filedir <- NULL
rm(filedir)

# Housekeeping
comp_dt <- NULL
panel_a_heading_input1a <- NULL
panel_a_heading_input1b <- NULL
panel_a_heading_input1c <- NULL
panel_a_heading_input1d <- NULL
rm(
    comp_dt,
    panel_a_heading_input1a,
    panel_a_heading_input1b,
    panel_a_heading_input1c,
    panel_a_heading_input1d
)

# APPENDIX G -------------------------------------------------------------------
# Here we calculate the numbers reported in Appendix G.

# IRS (2001): What is the MPE if the assumption is
# (i) delta = r = 5% (ii) delta = r = 2.5%?
a0 <- vector(mode = "numeric", length = 5)
a1 <- vector(mode = "numeric", length = 5)
a2 <- vector(mode = "numeric", length = 5)
for (t in 1:5) a0[t] <- lambda(r = 0.1, delta = 0.1, t = t, K = 50, T = 80) # average age in IRS (2001) is 50 #nolint
t <- NULL
rm(t)
for (t in 1:5) a1[t] <- lambda(r = 0.05, delta = 0.05, t = t, K = 50, T = 80)
t <- NULL
rm(t)
for (t in 1:5) a2[t] <- lambda(r = 0.025, delta = 0.025, t = t, K = 50, T = 80)
t <- NULL
rm(t)
alt_mpes_irs2001 <-
    c(round(-0.112 * mean(a0) / mean(a1), 3), round(-0.112 * mean(a0) / mean(a2), 3)) # nolint
rm(a0, a1, a2)

print(
    sprintf(
        "Under r = delta = 0.05, IRS (2001) MPE estimate becomes %s",
        alt_mpes_irs2001[[1]]
    )
)
print(
    sprintf(
        "Under r = delta = 0.025, IRS (2001) MPE estimate becomes %s",
        alt_mpes_irs2001[[2]]
    )
)

# CLNÖ (2017): What is the MPE if we assumed r = 2% and delta = 2.5% instead,
# using the simple back-on-the-envelope transformation of the existing estimate?
a0 <- vector(mode = "numeric", length = 80 - 49 + 1) # annuity factor when r == delta #nolint
a1 <- vector(mode = "numeric", length = 80 - 49 + 1) # annuity factor when r != delta #nolint
for (t in 0:(length(a0) - 1)) a0[t + 1] <- annuity_factor(r = 0.02, delta = 0.015, t = t, K = 49, T = 80) # average age in CLNÖ (2017) is 49 #nolint
t <- NULL
rm(t)
for (t in 0:(length(a1) - 1)) a1[t + 1] <- annuity_factor(r = 0.02, delta = 0.025, t = t, K = 49, T = 80) # nolint
t <- NULL
rm(t)
alt_mpe_clno2017 <- round(-0.284 * mean(a0[2:6]) / mean(a1[2:6]), 3)
rm(a0, a1)

print(
    sprintf(
        "Under back-of-the-envelope, CLNÖ (2017) MPE estimate is %s",
        alt_mpe_clno2017
    )
)

# Housekeeping
alt_mpes_irs2001 <- NULL
alt_mpe_clno2017 <- NULL
rm(alt_mpes_irs2001, alt_mpe_clno2017)
